library(Signac)
library(Seurat)
library(tidyverse)

# import peak mtx ----------
atac.norm <- Read10X_h5('~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/filtered_peak_bc_matrix.h5')
atac.norm[1:5,1:5]

atac.meta <- read_csv('~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/singlecell.csv') |>
  column_to_rownames('barcode')

chrom_assay <- CreateChromatinAssay(
  counts = atac.norm,
  fragments = '~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/fragments.tsv.gz',
  sep = c(":", "-"),
  min.cells = 10,
  min.features = 200
)

aobj <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = atac.meta
)

aobj[['peaks']]

granges(aobj)

# add genome annotation --------
## extract gene annotations from EnsDb packages
# BiocManager::install('EnsDb.Hsapiens.v86') # this is latest hg38 EnsDb
## create from gtf 
# annotations <- GetGRangesFromEnsDb(ensdb = ensdb.v112)

ensdb.hs109 <- EnsDb('Homo_sapiens.GRCh38.109.sqlite')

annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86)

seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg38"

# add the gene information to the object
Annotation(aobj) <- annotations

# QC ------
# compute nucleosome signal score per cell
aobj <- NucleosomeSignal(object = aobj)

# compute TSS enrichment score per cell
aobj <- TSSEnrichment(object = aobj, fast = F)

head(aobj@meta.data)

# compute blacklist frac via Signac 
aobj$blacklist_fraction <- FractionCountsInRegion(
  object = aobj, 
  assay = 'peaks',
  regions = blacklist_hg38_unified
)

# add blacklist ratio and fraction of reads in peaks
aobj$pct_reads_in_peaks <- aobj$peak_region_fragments / aobj$passed_filters * 100
aobj$blacklist_ratio <- aobj$blacklist_region_fragments / aobj$peak_region_fragments

DensityScatter(aobj, x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)

aobj$high.tss <- ifelse(aobj$TSS.enrichment > 3, 'High', 'Low')
TSSPlot(aobj, group.by = 'high.tss') + NoLegend()

aobj$nucleosome_group <- ifelse(aobj$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = aobj, group.by = 'nucleosome_group') +
  ggtitle('Nucleosome signal QC')

# qc violin plot
VlnPlot(
  object = aobj,
  features = c('nCount_peaks', 'TSS.enrichment', 'blacklist_fraction',
               'nucleosome_signal', 'pct_reads_in_peaks'),
  pt.size = 0,
  ncol = 5
)

# remove cells of bad quality
aobj <- subset(
  x = aobj,
  subset = nCount_peaks > 3000 &
    nCount_peaks < 30000 &
    pct_reads_in_peaks > 15 &
    blacklist_ratio < 0.05 &
    nucleosome_signal < 4 &
    TSS.enrichment > 3
)
aobj

# Normalization & dim reduc & clustering -------
aobj <- RunTFIDF(aobj)
aobj <- FindTopFeatures(aobj, min.cutoff = 'q0')
aobj <- RunSVD(aobj)

# LSI component 1 will strongly correlate with seq depth
# So we will exclude it from downstream analysis
DepthCor(aobj)

aobj <- RunUMAP(object = aobj, reduction = 'lsi', dims = 2:30)
aobj <- FindNeighbors(object = aobj, reduction = 'lsi', dims = 2:30)
aobj <- FindClusters(object = aobj, verbose = FALSE, algorithm = 3)
DimPlot(object = aobj, label = TRUE) + NoLegend()

# gene activity matrix ---------
# cost ~5 min
gene.activities <- GeneActivity(aobj)

# add the gene activity matrix to the Seurat object as a new assay and normalize it
aobj[['RNA']] <- CreateAssayObject(counts = gene.activities)
aobj <- NormalizeData(
  object = aobj,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(aobj$nCount_RNA)
)

DefaultAssay(aobj) <- 'RNA'

FeaturePlot(
  object = aobj,
  features = c('KRT1','TRPV3','KRT14','CD3D'),
  pt.size = 0.1,
  max.cutoff = 'q95'
)

skin.marker <- list(
  kera = c('KRT1','KRT10','KRTDAP','KRT14','KRT5','COL17A1'),
  fibroblast = c('FBN1','LUM','COL3A1'),
  macrophage = c('CD163','MRC1','MARCO'),
  T.cell = c('CD3E','CD8A','IL17A','FOXP3'),
  DC = c('CD1A','CD207')
)

aobj |>
  DotPlot(c(skin.marker,'TRPV3'), cluster.idents = T) +
  RotatedAxis()

aobj |>
  DotPlot(c(late.kc,key_cytokine,kegg_mva,'TRPV3'),
          cluster.idents = T,
          assay = 'RNA') +
  RotatedAxis()

# Find differentially accessible peaks ------
# change back to working with peaks instead of gene activities
DefaultAssay(aobj) <- 'peaks'

## simple, acceptable & quick
da_peaks.wilcox <- aobj |>
  FindAllMarkers(only.pos = T, latent.vars = 'nCount_peaks')

head(da_peaks.wilcox)

## use this param will cost >2h for one cluster
da_peaks <- aobj |>
  FindAllMarkers(test.use = 'LR', latent.vars = 'nCount_peaks',
                 only.pos = T)

head(da_peaks)

## identify closest gene from granges
marker.peak.cluster <- da_peaks.wilcox |>
  as_tibble() |>
  dplyr::filter(p_val_adj < .05) |>
  pull(gene)

marker.peak.close0 <- aobj |> ClosestFeature(marker.peak.cluster)

region.v3 <- marker.peak.close0 |>
  as_tibble() |>
  dplyr::filter(gene_name == 'TRPV3') |>
  pull(query_region)

da_peaks.wilcox |>
  as_tibble() |>
  dplyr::filter(gene == region.v3)

# plot genomic region -------
pbmc |> CoveragePlot(
  region = region.v3,
  extend.upstream = 0,
  extend.downstream = 80000,
  assay = 'peaks'
)
